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We show that the energy density spectrum of the primordial gravitational waves has characteristic 

features due to the successive changes in the relativistic degrees of freedom during the radiation 

era. These changes make the evolution of radiation energy density deviate from the conventional 

adiabatic evolution, pr oc a~*, and thus cause the expansion rate of the universe to change suddenly 

\ jr^ • at each transition which, in turn, modifies the spectrum of primordial gravitational waves. We 

f^ ' take into account all the particles in the Standard Model of elementary particles. In addition, free- 

(^ , streaming of neutrinos damps the amplitude of gravitational waves, leaving characteristic features in 

^SJ ' the energy density spectrum. Our calculations are solely based on the standard model of cosmology 

, , ' and particle physics, and therefore these features must exist. Our calculations significantly improve 

O . the previous ones which ignored these effects and predicted a smooth, featureless spectrum. 
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I. INTRODUCTION 

Detection of the stochastic background of primordial gravitational waves has profound imphcations for the physics 

of the earl y un iverse and the high energy physics that is not accessible by particle accelerators [HQjSSIMIaHBS 

^h ' llCt llll I12I |l3J . The basic reason why relic gravitational waves carry information about the very early universe is that 

f^ ' particles which decoupled from the primordial plasma at a certain time, t ~ idee, when the universe had a temperature 

^O of r ^ Tdec, memorize the physical state of the universe at and below Tdoc- Since gravitons decoupled below the 

^p Plank energy scale, the relic gravitons memorize all the expansion history of the universe after they decoupled and 

(^ thus would probe deeper into the very early universe. Gravitational waves act therefore as the time machine that 

rS 1' allows us to see through the entire history of the universe. Another example of relic species is the cosmic microwave 

1. background (CMB) photons, which decoupled from matter at T '^ 0.3 eV and can trace the physical state of the 

^ universe back to 0.3 eV. On the other hand, the primordial gravitational waves carry information on the state of the 

"^ much earlier universe than the CMB photons do. 

C^ The purpose of this paper is to study the evolution of primordial gravitational waves through changes in the physical 

L^ ' conditions in the universe within the Standard Model of elementary particles and beyond. For instance, the quark 

• i-H , gluon plasma (QGP) phase to hadron gas phase transition causes a sharp feature in the gravitational wave spectrum. 

^^ ' The change of the number of relativistic degrees of freedom affects the Hubble rate by reducing the growth rate of 

H I the Hubble radius during the transition. Thus, the rate at which modes re-enter the horizon is changed during the 

■ - - < transition and a step in the spectrum appears at frequencies on the order of the Hubble rate at the transition. For 

the QGP phase transition this frequency is ~ 10~^ Hz today and the correction is about 30% 14]. Other large 

drops in the number of relativistic degrees of freedom occur at electron positron annihilation and possibly at the 

supersymmetry (SUSY) breaking. Since the gravitational wave spectrum is sensitive to the number of relativistic 

degrees of freedom, one can search for evidence of supersymmetry in the very beginning of the universe by looking at 

the relevant frequency region (^ 10"'^ Hz). For energy scales lower than neutrino decoupling (~ 2 MeV [l^) we shall 

also account for the damping effect from neutrino free-streaming. 

The primordial gravitational wave spectrum will also provide us with information about inflation. The energy 
scale of inflation is directly related to the amplitude of the spectrum. The modes which re-entered the horizon 
during the radiation dominated epoch show a nearly scale invariant spectrum if we do not consider the change of 
the effective number of degrees of freedom. Typically the amplitude of the spectrum is of order 10^^^ for 10^^ GeV 
inflation ener gy scale in such a frequency region. Inflation ends when the inflaton decays into radiation and reheats 
the universe ^Q, [13, uM- The energy scale of reheating could be seen from the highest frequency end (~ fcrh) of 
a nearly scale invariant energy density spectrum of the primordial gravitational waves. The lowest frequency mode 
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observable today corresponds to the horizon size today, and the interval between the lowest frequency and fcrh would 
give the number of e-holdings, which tells us the duration of inflation between the end of inflation and the time at 
which fluctuations having the wavelength of the current horizon size left the horizon during inflation. The slope of 
the spectrum provides the power-law index of the tensor perturbation, ut |lCt Il2l|. tit = corresponds to a scale 
invariant power spectrum from de Sitter inflation. In a large class of inflationary models JtitI is not zero but much 
smaller than unity, and its determination constrains the inflationary models. As the effect of tit has been investigated 
by many authors, e.g. [13, lU H^ U^ 123 > ^'^'i i^ easy to include, we shall assume de Sitter inflation {nx = 0) 
throughout this paper. Our result is general and easily applicable to any kind of models which produce primordial 
tensor perturbations, (e.g. Ekpyrotic models J2l|'). 

The primordial gravitational waves not only test and probe the physics of inflation and reheating, but also can 
provide the tomography of the standard model of particle physics and models beyond. The study of its spectrum 
enables us to probe the very early universe in a truly transparent way. The goal of this paper is to show how the 
constituents in the early and very early universe would affect the primordial gravitational wave spectrum, which is 
observable in principle and may be observable in the future by the next generation observational projects, such as 
the Big Bang Observer (BBO) proposed to NASA [22| and the DECIGO proposed in Japan '23\. We present a new, 
rigorous computation of the primordial gravitational wave spectrum from de Sitter inflation with the Standard Model 
of particle physics. It is easy to extend our results to non-de Sitter (e.g., slow-roll) inflation models. 

The outline of this paper is as follows. In Sec. m basics about the primordial gravitational waves from inflation 
are reviewed. In Sec. IIIII a crucial quantity during radiation domination, the effective relativistic degrees of freedom, 
g■^,, is introduced and related to the primordial gravitational wave spectrum in heuristic and intuitive manners to 
illustrate the underlying physics. In Sec. IIVI we give an improved calculation of the primordial gravitational wave 
spectrum in the Standard Model, employing de Sitter inflation. Our final results are summarized in Figs. ^ and |S1 
In Appendix ^ we give useful formulae for the Bessel type functions. In Appendix ^ we give analytical solutions 
of gravitational waves in some limiting cases. We define energy density of gravitational waves in Appendix |0 The 
effect of neutrino free-streaming on the spectrum is formulated and explained in Appendix^ The numerical solution 
to the integro-differential equation is also presented. In Appendix ^ we give more detailed analytical accounts of 
numerical solutions of gravitational waves when the effective number of relativistic species changes. Units are chosen 
as c = ?i = fcs = 1 and y/SirG is retained. Indices A, /i, z^, . . . run from to 3, and i, j, k, . . . run from 1 to 3. 
Over-dots are used for derivatives with respect to time throughout the paper. Primes are mainly used for derivatives 
with respect to conformal time, but sometimes with respect to arguments we are focusing on. Barred quantities are 
unperturbed parts of variables. 

II. WAVE EQUATION, POWER SPECTRUM, AND ENERGY DENSITY 

In this section we define the power spectrum, A^(/c), and relative spectral energy density, flfi{k), of the gravitational 
wave background. We do this because some authors use different conventions in the literature. For tensor perturbations 
on an isotropic, uniform and flat background spacetime, the metric is given by 

ds'^ = a'^{T)[-dT'^ + {S^j +h,j)dx'dx^], (1) 

gt„^ = a^{T){T]f,^ +hf,„), (2) 
where 

Vfj^i' = diag(-l, 1,1,1), hoo = hoi^O, |/ijj| < 1. (3) 

Here and after we shall work in the transverse traceless (TT) gauge, which leaves only the tensor modes in pertur- 
bations, i.e. hijj = and h^i = 0. In the linear perturbation theory the TT metric fluctuations are gauge invariant 
^. We shall denote the two independent polarization states of the perturbation as A = -f , x and sometimes suppress 
them when causing no confusion. We decompose hij into plane waves with the comoving wave number, |k| = k, as 

--^ f d^k 
h,,{r,^) =J2j^-^h,{T;k)e^^%, (4) 

where e^- is the polarization tensor and A = -I-, x. The equation for the wave amplitude, h\{T;'k.) = /lA.ki is obtained 
by requiring the perturbed metric [Eq. jSJl] to satisfy the Einstein equation to 0{h). One finds that SGij = SirGSTij 



^ In classic references l24ll25l . hij = 2H'j'ij and Ilij = piTTij for tensor perturbations, which are automatically gauge-invariant. 



in the linear order ^g| yields 



--h.j.y" ^8ttGU,j, (5) 

where Ily {t, x) is the anisotropic part of the stress tensor, defined by writing the spatial part of the perturbed 
energy-momentum tensor as 

Ttj = pgij + a^IIy , (6) 

where p is pressure. For a perfect fluid Xly = 0, but this would not be true in general. In the cosmological context, 
the amplitude of gra vitational waves is affected by anisotropic stress when neutrinos are freely streaming (less than 
~ 10"'^°K) [23, 123, 123, yfl, IsJ, |33) ISS yJI • As we only deal with tensor perturbations, /i^, we may treat each component 
as a scalar quantity under general coordinate transformation, which means e.g. hij-^ — hij^^. The left-hand side of 
Eq. © becomes 

= ~'^'^ + i'^)'^'^ ~ ( T )^*^' ^''^ 

where Fq^^ — F^q = 0, F" — Sijda, and 5'-' = a^'^Sij have been used. Commas denote partial derivatives, while 
semicolons denote covariant derivatives in Eqs. ^ and Q. Transforming this equation into Fourier space, one 
obtains 

hx.M + (^) hx.M + (^) /^A,k = 167rGnA,k, (8) 

where the Laplacian V^ in the second term of |5j has been replaced by — fc^ in the third term of (jSJ. The second term 
represents the effect of the expansion of the universe. Using conformal time derivatives ( ' = ^), we may obtain 

h'U + (^) '^A.k + fc'/iA,k = IBTrGa^nA.k. (9) 

This is just the massless Klein-Gordon equation for a plane wave in an expanding space with a source term. Thus, 
each polarization state of the wave behaves as a massless, minimally coupled, real scalar field, with a normalization 
factor of y/lGirG relating the two. 

Next, let us consider the time evolution of the spectrum. After the fluctuations left the horizon, k <C aH, equa- 
tion would become 

''^^^^^ '' (10) 
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whose solution is 






(11) 



where A and B are integration constants. Ignoring the second term that is a decaying mode, one finds that /i^.k 
remains constant outside the horizon. Note that we have ignored the effect of anisotropic stress outside the horizon, 
as this term is usually given by causal mechanism which must vanish outside the horizon. Therefore, one may write 
a general solution of /iA,k at any time as 

ft,,k(r) ^ hll"'TiT,k), (12) 

where h^^ is the primordial gravitational wave mode that left the horizon during inflation. The transfer function, 
T{T,k), then describes the sub-horizon evolution of gravitational wave modes after the modes entered the horizon. 
The transfer function is normalized such that T{t, k) -^ 1 a.s k ^ 0. The power spectrum of gravitational waves, 
A\{k), may be defined as 

/r!k 
Y^l{r,k), (13) 



which imphes 

Al{r,k)^—Y,{\h^,^{r)f). (14) 

A 

Using equation (|12() , one may write the time evolution of the power spectrum as 

Al{T,k)^Al^„^[TiT,k)]\ (15) 



where 
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We have used the prediction for the amphtude of gravitational waves from de-Sitter inflation at the last equality, and 
i/jnf is the Hubble constant during inflation. One may easily extend this result to slow-roll inflation models. 

The energy density of gravitational waves is given by the 0-0 component of stress-energy tensor of gravitational 
waves: 

The relative spectral energy density, ri/i(T, fc), is then given by the Fourier transform of energy density, p/i(r) 
divided by the critical density of the universe, Pci{t) (see Appendix O for full derivation): 



dlnfc ' 



2 



Note that 0,,(A;) is often defined as flhiT,k) = i2^2j~fj^7^ k"^ [T{T,k)f = i2a^(T)g^(r) ^^('^' ^) ^^ ^^^ litera- 
ture [l3,[iy,l23- This definition is not compatible with the 0-0 component of stress energy tensor; however, it is a 
good approximation when the modes are deep inside the horizon, k ^ aH . Let us briefly explain a relation between 
these two definitions. The transfer function is usually given by Bessel type functions, T{x) — -^[Ajn{x) + i?y„(x)]. 
The conformal time derivative of the transfer function is thus given by -pT(x) = — ^[Ajn^i{x) + Byn+i{x)], 
where x = kr. Therefore, in the limit that the modes are deep inside the horizon, k 3> aH, one obtains 

^h{k) = ^2Slr)T^(r) [^'(^>^)]' ~ 12^1")^^ (r) ^'' [^i^^k)]^^ which agrces with the definition of nh{k) in [iSllllli. 
The difference between Jl/j and k^A\ would affect the prediction only at the largest scales, where both the overall 
amplitude and phases are different. (The phases are shifted by 7r/2.) 

Figure n shows a numerical calculation of equation H18|l for Q.,n = 1 — ^r, ^rh^ — 4.15 x 10^^, and h — 0.7. We 
ignored the contribution from dark energy, which is only important at the lowest frequency regime that we are not 
interested in in this paper. One may understand the basic features in this numerical result as follows. Energy density 
of gravitational waves evolves just like that of radiation inside the horizon, pft,(r, k) oc a""*, for k ^ aH . This implies 
that the relative spectral energy density, ^h{T, k), inside the horizon remains independent of time during the radiation 
era while it decreases as flhiT, k) (x a~^ during the matter era. Therefore, the modes that entered the horizon during 
the matter era later would decay less. As the low frequency modes represent the modes that entered the horizon at 
late times, fth{T, k) rises toward lower frequencies. On the other hand, flh{T, fc) at fc > 10^^^ Hz is independent of k. 
These are the modes that entered the horizon during the radiation era for which Qfi{T,k) was independent of time. 
After the matter-radiation equality all of these modes suffered the same amount of redshift, and thus the shape of 
^h{T,k) still remains scale-invariant at fc ^ 10^^^ Hz. 

These qualitative arguments may be made more quantitative by using the following analytical solutions of Q^ (t, k) 
for three different regimes (see Appendix IbI for derivation): 

f],(r<req,fc>fceq) = to^r" k'inikr)]', (19) 

17,(T>r,q,fc>fceq) = to'^Ta^ fc'^ [A{k)MkT) + B{k)y,{kT)f , (20) 

-L^-fJn ar\ T 



O ( -^ h ^ h \ — h,prim^ j^ 

'■l'h(T > Teq, K <. Koq) — ^ r. tt2 3 '*' 
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FIG. 1: The primordial gravitational wave spectrum at present, r — tq, as a function of the comoving wavenumber, k (or kc 
in units of Hertz). The frequency of gravitational waves observed today is related to k by /o — kc/2-K. The spectrum at large 
wavenumber is exactly scale-invariant as we have assumed de-Sitter inflation. In this figure we have not taken into account the 
effects of the change in effective relativistic degrees of freedom or neutrino free-streaming. 



where r^q is the conformal time at the matter-radiation equality, fee 
that entered the horizon at equahty, and Jq{x) — —ji{x) and 



31 (^) 



is the comoving wavenumber of the modes 
: —^^^ have been used to compute T'{t, k). 

(Spherical Bessel functions are given in Appendix 1X1) The first solution [Eq. (|19|l ] describes ri/i(T, k) during radiation 
era for the modes that entered the horizon before Teq. This solution is of course not relevant to what we observe today. 
(We do not live in the radiation era.) The second [Eq. H2UI) ] and third [Eq. H21|l ] solutions describe ri;j(r, fc) during 
matter era for the modes that entered the horizon before and after Tcq, respectively. The fc-dependent coefficients 
A{k) and B{k) are given in equation (|B9I) and (|B10|) . respectively. While the expression is slightly complicated, 
one can find that the second solution is independent of k when the oscillatory part is averaged out, which explains 
a scale- invariant spectrum at high frequencies, k > k^q ^ 10^^'' Hz. On the other hand, the third solution gives 
il;i(r, k) ex k~'^, which explains the low frequency spectrum. 

Figure m (and its extension to slow-roll inflation which yields a small tilt in the overall shape of the spectrum) has 
been widely referred to as the prediction from the standard model of cosmology. However, as we shall show in the 
subsequent sections, the standard model of cosmology actually yields much richer gravitational wave spectrum with 
more characteristic features in it. 



III. THE EFFECTIVE RELATIVISTIC DEGREES OF FREEDOM: 



It is often taken for granted that energy density of the universe evolves as p oc a~ during the radiation era. This 
is exactly what caused a scale invariant spectrum of i^h{k) at fc > kcq. However, p oc a^"' does not always hold even 
during the radiation era, as some particles would become non-relativistic before the others and stop contributing to 
the radiation energy density. 

During the radiation era many kinds of particles interacted with photons frequently so that they were in thermal 
equilibrium. In an adiabatic system, the entropy per unit comoving volume must be conserved p3a |: 



S{T) = s{T)a {T) = constant, 
where 

siT) = ^g„(r)T3. 



(22) 



The entropy density, s{T), is given by the energy density and pressure; s — {p + p)/T. The energy density and 



TABLE I: Particles in the Standard Model and their mass and helicity states 



particle 


rest mass [MeV] 


the number of helicity states: 


9i 


7 





2 




v,v 





6 




e+,e' 


0.51 


4 




fJ-'^,fi~ 


106 


4 




■K~^ ,TT^ 


135 


2 




tt" 


140 


1 




gluons 





16 




u,u 


5 


12 




d,d 


9 


12 




s,s 


115 


12 




c,c 


1.3x10-^ 


12 




T+,r~ 


1.8x10^ 


4 




bfi 


4.4x10^ 


12 




w+,w~ 


80x10^ 


6 




z 


91x10^ 


3 




H 


114x10^ 


1 




t,i 


174x10^ 


12 




SUSY particles 


~ 1 X lO** 


~ 110 





pressure in such a plasma-dominant universe are given by 



piT) 
P{T) 



—gJT)T\ 
IpiT), 



(23) 
(24) 



respectively, where we have defined the "effective number of relativistic degrees of freedom", g^, and g^,s, following |36l |. 
These quantities, g^, and g^,s, count the (effective) number of relativistic species contributing to the radiation energy 
density and entropy, respectively. One may call either (or both) of the two the effective number of relativistic degrees 
of freedom. Equation 122|) and H23|) immediately imply that energy density of the universe during the radiation era 
should evolve as 



-4/3 - 

p cx (7*5*s a 



(25) 



Therefore, unless g* and g^,s are independent of time, the evolution of p would deviate from p oc a~^. In other words, 
the evolution of p during the radiation era is sensitive to how many relativistic species the universe had at a given 
epoch. As the wave equation of gravitational waves contains {a' /a)h'^ j., the solution of h\,k would be affected by g^, 
and g^s via the Friedman equation: 



a'ir) 




-4/3 






Qr 



(26) 



Although the interaction rate among particles and antiparticles is assumed to be fast enough (compared with the 
expansion rate) to keep them in thermal equilibrium, the interaction is assumed to be weak enough for them to be 
treated as ideal gases. In the case of an ideal gas at temperature T, each particle species of a given mass, rrii — xiT, 
would contribute to g^ and g^,s the amount given by 



15 '■°° ^-2 

4 



TT 



(u2 _ a;f )i/2 



g*4T) = gt^ I '"" ^„ '^^ — u^du, 



, iT) _ 15 r- (.^-..f)V %^,, .^ 



[u^ - ^)du, 



(27) 
(28) 



1000 



100 



10 



1 

1e-06 0.0001 0.01 1 100 10000 1e+06 1e+08 

T [MeV] 

FIG. 2: Evolution of the effective number of relativistic degrees of freedom contributing to energy density, g*, as a function of 
temperature. The sohd and dashed lines represent <;« in the Standard Model and in the minimal extension of Standard Model, 
respectively. At the energy scales above ~ 1 TeV, gf^ = 106.75 and ^^ssm ^ 22O. At the energy scales below ~ 0.1 MeV, 
g, = 3.3626 and g,s = 3.9091; g« = g^s otherwise. 
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FIG. 3: Evolution of a' as a function of the conformal time. If g* and 3*3 were constant, p oc a * and a' would also be constant. 



where the sign is + for bosons and — for fermions. 

Here, gi is the number of heUcity states of the particle and antiparticle. Note that an integral variable is defined as 
u = E/T, where E = \/\p\^ + rm? ■ We assume that the chemical potential, /Xi, is negligible. One might also define a 
similar quantity for the number density. 



n{T) ^ ^g.„r3, 
where C(3) ~ 1.20206 is the Riemann zeta function of 3. Each species would contribute to 5*„ by 



g*n,i(T) = gi 



{u^-xlfl^ 



\du. 



(29) 



(30) 



■2C(3)7,, e«±l 

The effective number relativistic degrees of freedom is then given by the temperature-weighted sum of all particle 



contributions: 



9*(.T) = J29*AT) TfT 



T, 



>(r) = E' 



sAT) [^ 



9*n{T) =^g*nAT) TfT 



T, 



(31) 



where we have taken into account the possibihty that each species i may have a thermal distribution with a different 



T-v. Neutrinos are cooler 



L^ 



temperature from that of photons. The most famous example is neutrinos: T^ 

than photons at temperatures below MeV scale due to photon heating from electron-positron annihilation. 

Fig. 12 shows the evolution of g* as a function of temperature. We have included all the particles in the Standard 
Model of elementary particles, as listed in table |I| (Note that we assume that the mass of Higgs bosons is 114 GeV, 
which is the current lower bound from experiments.) We neglected hadrons whose mass is heavier than pions. In 
addition to the particles in the Standard Model, one may also include particles in supersymmetric models. Super- 
partners in the minimal extension of supersymmetric Standard Model (MSSM) would carry almost the same g^, as 
that carried by particles within the Standard Model. Fig. |31 shows the evolution of a'. If g^ and g^s were constant, 
a' would also be constant during the radiation era; however, the evolution of a' reveals a series of jumps due to the 
change in g^, and g^,s. 

Interactions between particles would change the ideal gas result obtained above, and one cannot use equation H27|l . 
(|28|l and (|31|l to calculate g* or g^s- Instead, one needs to extract g* and g^,s directly from energy density and entropy 
which would be calculated using detailed numerical simulations of particle interactions. For example, above the critical 
temperature of Quark Gluon Plasma (QGP) phase transition, most of g* is carried by color degrees of freedom. The 
dominant correction therefore comes from the colored sector of the Standard Model, whereas corrections from the 
weak charged sector are suppressed by the masses of weak gauge bosons. Since physics of QCD correction is still 
uncertain and beyond the scope of this paper, we shall ignore this effect and treat it as an ideal gas case. The effects 
of particle interactions on g^ have been investigated by 15^ .37 , ,2M ■ 



A. Heuristic argument based on background density 

Before presenting the full numerical results, let us briefly describe how g^ and g^s would affect the shape of Q^iTQ, k). 
In Sec. ^ we discussed how the expansion of the universe would affect ^.^{tq, k). While energy density of the universe 
during the radiation era is affected by 5* and g^,s as pcr oc g*^*^ a~*, energy density of gravitational waves always 
evolves as p;i(r, k) ex a"'* inside the horizon, k 3> aH , regardless of g* or g^s- (Gravitons are not in thermal equilibrium 
with other particles.) This difference in the evolution of ph and pcr significantly modifies a scale-invariant spectrum 
of Jlft(To, A;) at A: > Acq- 

Let us consider a gravitational wave mode with k which entered the horizon at a given time. The < '''cq and 
temperature, T = The, during the radiation era. After the mode entered the horizon the amplitude of this mode 
would be suppressed by the cosmological redshift. The relative spectral density at present would then be given by 



ri/i(ro,fc > Acq) = ri/i(rhc,fc)f7ro 



ff*s(T'hc) 



g*so 



-4/3 



g*iThc) 



g*o 



(32) 



where fir denotes the relative energy density of radiation and the subscript "0" denotes the present-day value. This 
equation helps us understand how 5* and g^,s would affect J7/i(To,fc). For a given wavenumber, fc, there would be 
one horizon-crossing epoch, The- The amount by which the relative spectral energy density of that mode would 
be suppressed depends on g^, and g^,s at The. The mode that entered the horizon earlier should experience larger 
suppression, as g* and g^s would be larger than those for the mode that entered the horizon later. (The effective 
number of relativistic degrees of freedom is larger at earlier times — see Figure |21) As g^ and g^s are equal for 
T > 0.1 MeV and nearly the same otherwise (g* = 3.3626 and g*s = 3.9091 for T < 0.1 MeV), we expect that 
suppression factor is given by ((7*/5*o)~^ to a good approximation. The modes that entered the horizon during the 
matter era should not be affected by g* or g^s, as they do not change during the matter era. 



B. More rigorous argument using analytical solutions 



In this subsection we derive equation (|32|l using a more rigorous approach. Let us go back to the wave equation 
[Eq. (PI], and rewrite it using a new field variable, fik = ahk- 

fi'k + (fc2 - — Lfc = mnGa^Ilk. (33) 



Note that we have suppressed the subscript for polarization, A. To find a solution for /i^, we must solve the Friedman 
equation as well: 

4)' = ^P^--^o^^(^)"'7-)" (34) 

, / X 1/2/ X -2/3/ X -2 

a V.9*o/ \g*so/ \ao, 

pT / \ 1/2 / \ —2/3 

ao J TO \9*oJ \g*o 

where the subscript "0" denotes some reference epoch during the radiation era. (While "0" means the present epoch 
in the other sections, we use it to mean some epoch during the radiation era in this section only.) To proceed further, 
we need to specify the evolution of g^ and g.^,s. While we have numerical data for the evolution of these quantities, 
we make an approximation here to make the problem analytically solvable. Since (7*(r) decreases monotonically as 
the universe expands, one may try a reasonable ansatz, g* oc r"^", to obtain analytical solutions. We shall also 
assume g* = g^,s and Hk — for simplicity in this section. (At temperatures below 2 MeV, free-streaming of neutrinos 
generates anisotropic stress, 11^; ^ 0. Also, the temperature of neutrinos is different from that of photons below 
electron-positron annihilation temperature, and thus g* ^ g^s below ^ 0.51 MeV.) This model gives 

a" 9*''^^ I 5* I " ,-, , X ,or\ 

— ~ -7-~— -(l + n), (35) 

where the primes denote derivatives with respect to r, and the term -t- (1 + n) has been neglected in the last line, 
assuming r <C tq. This form of a" / a allows us to find an analytical solution to equation (|33|l : 

where A{k) and B{k) are the normalization constants that should be determined by the appropriate boundary con- 
ditions. Note that n — and n — 1 correspond to the solutions for the radiation era and the matter era, respectively. 
Let us consider a model of the radiation-dominated universe in which there was a brief period of time during which 
g* suddenly decreased as a power-law in time, 17* oc t~^". Outside of this period g* is a constant. Suppose that g^ 
changed between t = T2 and ri > T2. (The change in g^, began at r = T2 and completed at ti.) The modes that 
entered the horizon after ri do not know anything about the change in g^,. The solution for such modes is therefore 
given by the usual solution during the radiation era, 

h^ir) = hl""^joikr). (37) 

How about the modes that entered the horizon before r2? The solution for such modes is given by 

hi^ir) = hr'-joikr) (r < r,), (38) 

J prim 
K{r) = ^[A{k)jn{kT) + B{k)yn{kT)] (r2<T<Tl), (39) 



4"(r) = hl""^ (^aJ [C{k)j,{kT) + D{k)y,{kT)] (n < r). 



(40) 



It is convenient to define T.^, = (ti + T2)/2 and At = ti — T2 to characterize the time of transition and its duration, 
respectively. Here, the superscript "in" denotes the modes that have already been inside the horizon at r*, while 
"out" denotes the modes that are still outside the horizon at t*. The coefficients, A{k), B{k), C{k), and D{k), are 
given by Eqs. (jElll - (|E4p in Appendix El By taking a ratio of equation H40|l and H37f) . we can find the amount of 
suppression in h^^{T > ri) relative to h'^^{T > ti): 

K{r>T^) _ (t2\ ^c{k)+D{k)y^{kT)ljo{kT)]^(^-^\ [C{k) + D{k)l (41) 



where we have ignored the oscillatory part oi yQ^kr) / JQ{kT) . While C(fc) and D{k) have fairly cumbersome expressions, 
the sum of the two has a simple limit, [C{k) + D{k)]'^ -^ 1, for At -^ 0, regardless of the value of n (see Appendix IE|| . 
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The energy density in gravitational waves then reflects the effect from the change of g* as 



»k"(r>Ti,fc) 

r!™t(r>Ti,fc) 



Kir > n) 



/ir'(T" > n) 



At 
l-2n , (42) 



where we have used the sub- horizon limit for il/i(fc) and At <^T2 < t* < ti. On the other hand, g^ cc t gives 

ffr* 9*{ri) \t,+At/2J n ^ ' 

Hence, combining Eas. (|42|l and ()43|l . we finally obtain the desired result 



'1/3 

f^r*(fc) V.9r* 



^r(fc) 



(44) 



for At <S^ r*. This result agrees with equation (|32|l . which was obtained in the previous section fSec. IIIlX)l using a 
more heuristic argument. (Note that we have assumed 5* = 5*^ in this section). In Eq. (|32|) there is an extra factor 
ilr-Oi which represents the time evolution of flh from matter-radiation equality to the present epoch. We do not have 
this factor in equation (|44|l . as both fi™ and fi™* are evaluated during the radiation era. 

IV. PREDICTION FOR ENERGY DENSITY OF GRAVITATIONAL WAVES FROM THE STANDARD 

MODEL AND BEYOND 

In Sec. lIIII we have described how the evolution of the effective number of relativistic degrees of freedom would affect 
the shape of relative spectral energy density of primordial gravitational waves at present, ^^(to, k). In this section we 
present the full calculation of Qh^To, fc), numerically integrating the wave equation together with the numerical data 
of g* and g^,s (see Figure EJ. 

Before we do this, there is another effect that one must take into account. While we have ignored anisotropic stress 
on the right hand side of the wave equation so far, free-streaming of relativistic neutrinos which have decoupled 
from thermal equilibrium at T 5^, 2 MeV significantly contributes to anisotropic stress, damping the amplitude of 
primordial gravitational waves |3ll l32l | . Calculations given in Appendix ^ show that neutrino anisotropic stress 
damps flh{To, k) by 35.5% in the frequency region between ~ 10""'^^ and ~ 2 x 10~^" Hz. The damping effect is much 
less significant below 10"^^ Hz, as this frequency region probes the universe that is dominated by matter. One may 
understand this by looking at the right hand side of Eq. IID23I) . Anisotropic stress is proportional to the fraction of 
the total energy density in neutrinos, fi>{T), which is very small when the universe is matter dominated. 

We show the results of full numerical integration in Figure 0] and [S] The latter figure is just a zoom-up of interesting 
features in the former one. We find that 17;i(to, k) oscillates very rapidly as sin (fcr -I- (p), where ip is a. phase constant. 
The cross term, sin fcr cos fcr, appeared as a beat in Fig.^ while they are too small to see in Fig^ From observational 
point of view these oscillations will not be detectable, as observations are only sensitive to the average power over a 
few decades in frequency. 

The damping effect due to neutrino free-streaming is evident below 2 x 10^ ^*', while one might also notice a minor 
wiggly feature at around 5 x 10^^° Hz. This feature is actually artificial. We implicitly assumed an instantaneous 
decoupling of neutrinos from the thermal plasma at T^dec — 2 MeV, which resulted in the surface of decoupling 
that is extremely thin. This gave rise to dips and peaks corresponding to the waveform of gravitational waves at the 
decoupling time. (The envelope shape is somewhat similar to — ji(fc) at around 5 x 10"^" Hz; more details are given in 
Appendix^) Physically speaking, however, the last scattering surface of neutrinos is very thick, unlike for photons. 
(There is no "recombination" for neutrinos.) Therefore, the oscillatory feature would be smeared out when thickness 
of the decoupling surface is explicitly taken into account. To do this, one would need to solve the Boltzmann equation 
for neutrinos separately, including the effect of neutrino decoupling. 

The effect of evolution of 5* and g^,s is also quite prominent. For example, big changes in g^, would occur at 
the electron-positron annihilation epoch, ~ 0.51 MeV (~ 2 x 10~^^ Hz), as well as at the QGP to hadron gas phase 
transition epoch, ^ 180 MeV (~ 10~^ Hz) within the Standard Model. The gravitational wave spectrum is suppressed 
by roughly 20% and 30% above the electron-positron annihilation and QGP phase transition scale, respectively. If 
supersymmetry existed above a certain energy scale, e.g., ~ 1 TeV (~ 1 x 10~'*Hz), the spectrum would be suppressed 
by at least ^ 20% (for N=l supersymmetry) above that frequency. We also find additional features at the QGP 
phase transition scale, ^ 10^^ Hz, similar to the features at ~ 5 x 10~^*^ Hz caused by our assumption about 
instantaneous decoupling of neutrinos. The feature at the QGP phase transition is nevertheless not artificial — as 
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FIG. 4: The primordial gravitational wave spectrum at present, flh{To,k)/10~^^ , as a function of the comoving wavenumber, k 
(or kc in units of Hertz). The frequency of gravitational waves observed today is related to fe by /o = kc/2-K. We have assumed 
a scale-invariant primordial spectrum and Qrn. ~ 1 — f2r, ^r = 4.15 x 10~^h~^, h = 0.7, and Ei-nt ~ 10^'' GeV. We have included 
the effects of the effective number of relativistic degrees of freedom and neutrino free-streaming. The dashed line shows the 
envelope of the previous calculations which ignored the change in the number of relativistic degrees of freedom and neutrino 
free-streaming (Fig.0. 



the QGP phase transition is expected to have happened in a short time period, the instantaneous transition would 
be a good approximation, unhke for neutrinos. 

One may approximately relate the horizon crossing temperature of the universe to the frequency of the gravitational 
waves [13, |41|. The horizon crossing mode, /chc — aiic-ffhc, is related to the temperature at that time by H^^ = 



Stt^G 



*?*,hc 






Then using entropy conservation, 5*s,hcahc^hc ~ 9* 



sOOi 



90 y*,nc^hc 

from the temperature of the universe to the frequency of gravitational waves observed today 



oTg , one obtains the following conversion factor 



/o = 1.65 X 27r X 10 



^7 I The A 
IGeV; 



g*s(rhc) 

100 



-1/3 



9* (The) 



100 



1/2 



Hz, 



(45) 



which was derived in |l3jlil|- (If we take e = ^^ in |41|, their equation (156) agrees with the one above.) Throughout 
this paper we have been using the comoving wavenumber, k (or kc in units of Hertz), which is related to the 
conventional frequency by 27r/o = kc/ag, where uq is the present-day scale factor and c is the speed of light. We use 
k in this paper, rather than /o, as k is what enters into the wave equation that we solve numerically. 



V. DISCUSSION AND CONCLUSION 



We have calculated the primordial gravitational wave spectrum, fully taking into account the evolution of the 
effective relativistic degrees of freedom and neutrino free-streaming, which were ignored in the previous calculations. 
The formalism and results given in this paper are based on solid physics and can be extended to primordial gravitational 
waves produced in any inflationary models and high energy particle physics models. As is seen in Figs. 0] and El the 
spectrum is no longer scale invariant, but has complex features in it. Whatever physics during inflation is, one must 
include the evolution of the effective relativistic degrees of freedom and neutrino free-streaming. 

14] studied the gravitational wave spectrum at the QGP phase transition assuming the first order instantaneous 
model as well as the second order cross-over model, and found 30% suppression of the energy density spectrum, which 
is consistent with our calculation. [33 studied the efi^ect of entropy production from e.g., decay of massive particles in 
the early universe on the energy density spectrum. We have not included this effect in our calculations, as the late-time 
entropy production is not predicted within the Standard Model, [iy] studied the effect of changes in the equation 
of state of the universe on the energy density spectrum. While they included the effect of neutrino free-streaming. 
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FIG. 5: A blow-up of Fig. 2] Note that density of vertical lines shows density of sampling points at which we evaluate f2h(ro, k). 
The dashed line shows the envelope of the spectrum in the Standard Model of elementary particles. 



they did not include the evolution of g*. Instead, they explored general possibilities that the equation of state might 
be modified by trace anomaly or interactions among particles. They also considered damping of gravitational waves 
due to anisotropic stress of some hypothetical particles. Our calculations are different from theirs, as we took into 
account explicitly all the particles in the Standard Model and the minimal extension of the Standard Model, but did 
not include any exotic physics beyond that. 

Let us mention a few points that would merit further studies. At the energy scales where supersymmetry is unbroken 
(if it exists), say TeV scales and above, the number of relativistic degrees of freedom, g*, should be at least doubled, 
and would cause suppression of the primordial gravitational waves (Fig.|51for A^ = 1 supersymmetry). If, for instance, 
A^ = 8 is the number of internal supersymmetric charges, ~ 250 copies of standard model particles would appear 
in this theory. This would suppress the spectrum by 85% at the high frequency region (above ~ 10~* Hz) compare 
to the Standard Model, though the details depend on models. Since we still do not have much idea about a true 
supersymmetric model and its particle rest mass, the search for the primordial gravitational waves would help to 
constrain the effective number of relativistic degrees of freedom 5* (T) above the TeV scales. 

In an extremely high frequency region, fcj-h, the gravitational wave spectrum should provide us with unique infor- 
mation about the reheating of the universe after inflation. If the inflaton potential during reheating is monomial, 
V{(p) ex 0", the equation of state during reheating is given by p^ — a{n)p^, where a{n) = ■^^^. Since the equation 
of state determines the expansion law of that epoch, one obtains the frequency dependence of the gravitational wave 
spectrum as flh oc fc("~'')/("~i). In an extremely low frequency region (below ~ 10~^^ Hz), on the other hand, dark 
energy dominates the universe and affects the spectrum. Acceleration of the universe reduces the amplitude of gravi- 
tational waves that enter the horizon at this epoch; however, we will not be able to observe modes as big as the size 
of the horizon today. 

The signatures of the primordial gravitational waves may be detected only by the CMB polarization in the low 
frequency region, ^ 10~^^ Hz. For the higher frequency region, however, direct detection of the gravitational waves 
would be necessary, and it should allow us to search for a particular cosmological event by arranging an appropriate 
instrument, as the events during the radiation era are imprinted on the spectrum of the primordial gravitational 
waves. 



APPENDIX A: SPHERICAL BESSEL TYPE FUNCTIONS 



We present some formulae for Bessel type functions used in this paper. 



dx 



Znix) 



Zn+lix) 



— [x"+hnix)] =.t"+1z„_i(x). 



(Al) 
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where Zn{x) can be spherical Bessel functions, spherical Neumann functions, Bessel functions, and Neumann functions. 
Spherical Bessel functions and spherical Neumann functions are related by 



ynix) = i-ir+'j-n-lix) 



Their asymptotic forms are 



jn{x) 



sin [x — m:/2) 



Vn{x) 



cos [x — mi/2) 



(A2) 



(A3) 



for a; ^ 1. If n is even, j,i(a;) « ±Jq(x) and yn{x) ~ ±yo(x). If n is odd, jn{x) ~ ±ya{x) and j/„(a:) ~ ±Jq{x). The 
first and second kinds of spherical Hankel functions are defined as 



h^nHx) = jn{x) + iy„(x), h^n\x) = Jn{x) - iy„(x) 

Using elementary functions, we have 

sin x 
jo[x) = — , 

ji{x) = 



sma; 



i2(a;) = - 
x 

yo{x) = 
yi{x) = 



X 



3 A . 3 ■ 

-^ — 1 ] sm X cos X 

X'^ I X 



cos a; 



2/2(2;) = 

X 



— cos X + sm X 

X 



3 A 3 . 

— — 1 I cos a; H — sm a; 



X \ X ' 






a: \ a; ' 



(A4) 

(A5) 
(A6) 

(A7) 
(A8) 
(A9) 

(AlO) 

(All) 

(A12) 



APPENDIX B: ANALYTICAL SOLUTIONS OF WAVE EQUATION 

In this Appendix we shall discuss solutions of the equation of motion [Eq. Q]. While we assume Ily = in this 
Appendix, we shall treat By 7^ in Appendix IdI Imposing appropriate boundary conditions |42l |. one obtains simple 
analytical solutions for tensor modes of fluctuations in the inflationary (de Sitter), radiation dominated (RD) and 
matter dominated (MD) universe, as 



hk{T) 



fi 



v'167rG / I, 1 _ji 

--VSTrGkh^j^\kT)a(k) 



hUr) - [jo{kT)]hP^ 



kr 



iprim 



inflation, 
RD, 

MD, 



(Bl) 
(B2) 
(B3) 



where a(k) is a stochastic variable satisfying (Q;(k)a*(k')) = (5''(k — k'), and spherical Bessel-type functions are given 
in Appendix 1X1 We classify wave modes by their horizon crossing time. The; 



, I 7 J > ^oq the modes that entered the horizon during RD: The < Tcq 
1 < fcoq the modes that entered the horizon during MD: The > ''cq 



(B4) 
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where r^q denotes the time at the matter-radiation equality, and The denotes the time when fluctuation modes crossed 
the horizon, /cThc = 1- Notice that [^.^(t)^ for each solution (|B1I) - ljB3|) does not depend on time (= \h'^""^\'^) at the 
super- horizon scale, |fcr| ^ 1. 

The tensor mode fluctuations from the inflationary universe left the horizon and froze out. Its dimensionless 
spectrum is given from Eq. (|B1(I as 






TT \ mpi J 1'k'^ 

where i/inf is the Hubble parameter during inflation and t = — l/(aiJinf) is used in the second equality. Note that 
the conventional factor 4 is from / ^^\{k) = {h.jh'^) = 2[{\h+\'^) + {\hy,\^)] = A\h\^, where \h+^k\ ^ l^x.fcl = \h\ is 
assumed [4^. From the Friedman equation during inflation, one obtains H?^f f« ^^^ ]/(</>), which gives A^^^^^ « 

lOy(0)/mpi; thus Aj^^^^^ is sensitive to the shape of inflaton potential [iflli3- The dimensionless spectrum HB5|) 
is nearly independent k. This is the famous prediction of the inflationary scenario known as a nearly scale invariant 
spectrum. As long as we consider de Sitter inflation, the spectrum is exactly scale invariant, i.e. oc fc" as </> is at rest. 
Using the transfer function [Eq. (|12|) ]. we obtain the time evolution of the amplitude of gravitational waves as 

T{t < Toq, k > fcoq) = joikr), (B6) 

r(T > Teq, fc > fceq) = ^ [A(fc) ji (fcr) + B(%i (fcr)] , (B7) 

T 

T{T,k<k,q) = ^^^, (B8) 



where 



^(;.) ^ _^ cos2fcrcq ^ sin2fcTeq ^ ^^^^ 

ZKTcq ZKTcq [r^'^cq) 



1 cos2kTcq sin2A;T( 

(fc^ " (fcTeq)2 2fc^ 



eq 



Bik) = -l + TTZ^^-^^^^^^^-^^TTTT^- (BIO) 

Their conformal time derivatives are given as 

T'(t< Teq, fc>A:cq) = ~kji{kT), (Bll) 

r'(T> Teq, fc> fceq) = - ^ [A(fc)j2 (^t) + i?(%2 (fc^)] , (B12) 

r'(T,fc<fceq) = -MaM. (B13) 

T 

Eqs. (|B6(I and (|B7|) are the evolution of modes which entered the horizon during the radiation era, while Eq. I)B8|I 
is the evolution of modes which entered the horizon during the matter era. Coefficients A(k) and B{k) are obtained 
by equating a solution (|B6|) with (|B7|I and their first derivatives [ (|B11|) and (|B12p ] at the matter-radiation equality. 
The transfer function for the intermediate regime, Eq. (|B7|I . can be calculated numerically so that the two other 
limiting solutions match smoothly (See Fig.^. If the wavelength of the gravitational waves is much shorter than the 
duration of the cosmological transition, a WKB approximation may be appropriate |32l |44|. Here we just assumed 
the instantaneous transition to illustrate the main point. The analytical solutions as well as numerical solutions are 
presented and compared in Fig. El and [3 The higher fc-modes enter the horizon earlier, and their amplitudes are 
damped more by the cosmological redshift. 

APPENDIX C: THE RELATIVE SPECTRAL DENSITY: Qnik) 

In this Appendix we shall define the energy momentum tensor of gravitational waves following the argument and 
the definition in §35.7 and §35.13 of ^45]. The Ricci tensor for the metric of the form given in Eq. (^ may be expanded 
in metric perturbations, h: 

R^.u = i?,,. + Ri'} + i?^j + 0{h'), (CI) 
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FIG. 6: Numerical solutions of tensor perturbations. The solid, dashed, and short-dashed lines show the high, medium, and 
low frequency modes, respectively. The higher fe-modes enter the horizon earlier, and are damped more by the cosmological 
redshift. Vertical lines define the horizon crossing time for each fc-mode. 
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FIG. 7: Comparison between numerical solutions and analytical solutions of tensor perturbations. The dashed and short-dashed 
lines show numerical solutions of the high and low frequency modes, respectively. The higher fc-modes enter the horizon earlier, 
and thus the numerical solution is well approximated by the analytical solution during the radiation era, xi^''') ~ joikr) (solid 
line). On the other hand, the lower fc-modes enter the horizon much later, and thus the numerical solution is close to the 
analytical solution during the matter era, xi^T) ~ 3ji{kT)/kT (dotted line). 



where R 



(1) 



0{h) and b!-^^ 



'fiiy 



0{h^ 



For the vacuum field equation, R^^, = 0. As the Einstein equation is non-linear, R^^ is in general not linear in h^^,. 
The linear term in Eq. (JC1|) must obey the vacuum equation. 



i?(i) = 



(C2) 

This is an equation for the propagation of the gravitational waves, which corresponds to Eq. @ or more generally to 
Eq. (|D23|I in the FRW universe. The remaining part of Rf^^, may be divided into a smooth part which varies only on 
scales larger than some coarse-graining scales, 



Rf.. + (41^) 



0, 



(C3) 
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and a fluctuating part which varies on smaller scales 

up to the second order in h^^. Here, i?J,J"°" '"''''' is defined by Eq. IJC4|I and represents the nonlinear correction to 
the propagation of h^j^^, Eq. IJC2|I . which gives h^^i, -^ h^^ + j^v, where j^y ^ 0{h'^) A5\. Eq. IJC3|) represents how the 
stress energy in the gravitational waves creates the background curvature. The Einstein equation in vacuum is then 

G^. = Rf.. - \r9^.u = 8^GT^^^), (C5) 

r(«^) ^ -^ ((4^)) - \g,AR''')) , (C6) 

where T^^, is a definition of the energy momentum tensor for the gravitational waves and ( ) denotes an average 
over several wavelengths. The importance of the effective energy momentum tensor is that it tells us how backreaction 
from energy density of gravitational waves would affect the expansion law of the background universe. Note that the 
effective energy momentum tensor defined by Eq. I|C6|I is different from that defined by the Neother current of the 
Lagrangian density, r^jsother ^ _2 ss^^ ^ where S'^^^ is the second order perturbation in the Einstein-Hilbert action. 
These definitions coincide only deep inside the horizon. Note also that in the notation of |45|. G = 1, but in our 
notation, h = c = 1, G = rripj'^ , where mp\ is the Plank mass. Since (i?*^^-') = p5| . 

where | is the covariant derivative with respective to background metric, g^^. Note that we have employed the 
transverse-traceless (TT) gauge. In linear theory we neglect higher order terms in the energy momentum tensor. 
The energy density of gravitational waves, ph, is defined by the 0-0 component of the energy momentum tensor. 



p,^(r)^T^^''> = ^^-^{h,,h^^), (C8) 

where /ly is in the TT gauge. There are only two independent modes for gravitational waves; 

f h+ /ix \ 

(C9) 

where + and x denote two independent polarization modes and their propagation direction is taken in z direction. 
Hence, 




^ :{h'l + h'l: 



167rGa2 

1 f d\ [ (Pk' 



167rGa2 J (27r)3 J (27r) 



((/i'+,k/i'+,k' +/'/x,k/i'x,k')e*'''+'''''"), (CIO) 



where Fourier transformation was done and /i*A,k = ^A,-k in the last step. For stochastic modes, the spatial average 
over several wavelengths, ( ), is equivalent to the ensemble average in k space; 

{h'xMh'x',^') = {2T:f5x,x'5^''\\^ + ^')\h'x.k\\ (Cll) 

where A = +, x. Using 1C10|I and IClip . we obtain 

1 /■ d^k 

It is reasonable to assume that the primordial gravitational waves are unpolarized, i.e. j/i+^fcp — |ft.x,feP- When- 
ever we express the time evolution of some quantities, it is convenient to express them in terms of the transfer 
function, T(fcr), and the primordial amplitude,A? ^j^, defined as ((T^ : 

Ph{T) = 32^ fd\nkAl^„^[rikT)f, (C13) 
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where 

k^ ^P, / FT- c'^^ 

2 A '^ \ipr'irn\2 ( -'^"ii 



Here, |ft,^'^*'"p is the ampUtude of gravitational waves outside the horizon, | kr |<C 1, during inflation. WeU inside 
the horizon averaging over several periods, the leading term of [T'(fcr)]2 is proportional to r~^ ex a~^ during the 
radiation era and ex r"'* oc a~^ during the matter era. Thus ph oc a^**, which is consistent with the fact that graviton 
is massless and thus relativistic. 

It is common to define the relative spectral density as the normalized energy density per logarithmic scale. 

n,ir,k) ^ P^, (C15) 

Ph[r,k) 



d\nk 



where Pcr{T) is critical density of the universe, and ph{T,k) denotes energy density of the gravitational waves per 
logarithmic scale. Inserting (|C13(I into (|C15|I . we obtain 

^^i^^k) = ^^'7" A r'{r,k)f- (C16) 

Recalling Friedman equation, H^ = SirGpc/S, ljC16|) becomes 

In this paper, we shall evaluate this quantity exactly within the Standard Model of elementary particles. For an 
analytical model, T'{t, k) is given by Eqs. IJB11|) - ljB13|l . 

APPENDIX D: COLLISIONLESS DAMPING DUE TO NEUTRINO FREE-STREAMING 

In this Appendix, we review the effect of collisionless particles on gravitational waves. Treating relativistic neutrino 
gas by classical kinetic theory, the linearized Einstein-Boltzmann equation Q can be written as an integro-diffcrcntial 
equation l|D23|) . The derivation of this integro-differential equation is given in the literature, for instance 27 , 29 , 30, 3J 
for both scalar and tensor modes, [23 for scalar modes, and will be reviewed briefly in this Appendix. 

At the temperature of '--^ 2 MeV, where neutrinos decoupled and became out of equilibrium with photons, electrons, 
or positrons, the number of effective relativistic species is g*{^ 2MeV) — 10.75.^ The free-streaming neutrino gas 
after their decoupling satisfies the collisionless Boltzmann equation, i.e. the Vlasov equation, 

^^^=0, (Dl) 

dt ^ ' 

where F{x, P) — F{P) + 5F{x, P) is a distribution function. The distribution function of relativistic neutrinos is 
given by 

where gi, denotes the number of helicity states for neutrinos and anti-neutrinos. Here, P^ = ^j^ and P" — ^/gtjP^Pi , 
which is implied by the constraint for relativistic particles; 

ff^.P^P" = 0. (D3) 



We have assumed instantaneous decoupling of neutrinos, but this is not true in general. 
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Therefore, there are only three independent components of the momentum vector. One can also relate P' with 
Pn = -P° as 



P^ = ± 



7'Po 



1 



hjkl^-i^ 



(D4) 



where 7* = 7i's are directional cosines and P'^ is the energy of neutrinos. We chose positive sign convention for 
P° = ^. Note that S^jYj' = 1, and P* = C7'Po, where the coefficient, C, is obtained from Eq. HD3|) : 



PnP° 



a^pipi +a^hijP'P\ 



= ^oJ 

= -(Po)2 + a'c2p2 + a2/iy7VC'^o, 



1 = 0^^2(1 + /i,^.yy). 

We consider tensor perturbations. Eq. IJD1|) can be expressed as 

dF{t,x\-f\P°) _ dF dx' dF dP" dF d-y' dF 
dt ~ 'dt^dt dx' "^ dt 9P0 ^ "dT^y 



= 0. 



(D5) 



The last term is negligible in the linear perturbation theory, as j~ is of the first order in perturbations and 7* 



-hh3k,^l^l'"■ 



d'y^ 



For the second term ^ is of the first order in perturbations and 



Using Eq. (|D4|1 . one obtains 



in the leading order, as F does not depend on x'; thus, Mr is a perturbation. 
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For the third term we use the geodesic equation 

rfP^ 



dX 



T^af^ 



r^„flP"P'3, 



9' 



.liv 



dxP dx°' dx'' 



(D8) 
(D9) 



The time component of the geodesic equation is 

dtdP^ 



dX dt 



-V\pP''pP, 



, dgau dga 



dxP dx" 



papp 



0\2 



-in 



2 dt 



(DIO) 



where (700 = — 1, goi — were used from the second line to the last line. Up to the first order in perturbations 

1 dpo a IdKj , , 

^7 7^ 



po dt 



a 2 dt 



(Dll) 



where we have used Eq. IJD4|) and neglected higher order terms. This equation describes the change in the neutrino 
energy as it propagates in a FRW universe with gravitational waves. The first term accounts for the redshift of 
energy due to an isotropic expansion. The second term tells us that neutrinos lose energy if —g^ > 0, or gain energy if 

-g^ < from gravitational waves. This energy flow from neutrinos to gravitational waves causes collisionless damping 
(Figs. IHlandO and amplification fFig. llO|l of gravitational waves. 



19 
Finally, by combining Eqs. (jDSp . ljD7p . and (jDlip . the Vlasov equation for the first order perturbations is obtained 

dF\ _d6F^^^dSF pod6Fa^ d^ldh^ 



as 

dt J nrstor^ior dt a dx^ dP° a dP"2 dt 

where F = F + 5F{t,x^ ,^^ ,P^) and 5F is a tensor type perturbation in a distribution function of neutrinos. The 
zeroth order Vlasov equation merely gives cosmological redshift, P" ex a^^, as explained above. Defining ji = Y^i/k 
and Fourier transforming Eq. IJD12I) . the first order Vlasov equation in the momentum space is given as 



where we have used 





dfk 
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a po dfk 
a dP° 


h^ 


j(i,x) 


A=+,x ■ 




SF 


= >; 



a dP° 2 dt ' 



(27r)3 

d^k 

(2^ 



A,fe(i,P°,/i)f7^g^,(x). (D15) 



Here, tensor harmonics Qfj{^) are solutions of the tensor Helmholtz equation; Q^..J'^{x.) + k^Q^j{x) = 0, diQ^j = 
ikiQ^j. They are symmetric, traceless, and divergenceless; Q^^ = Qj^, j^-'Qij = Qij^-^ = 0, where 7*^ = a^g*-' and | 
denotes the covariant derivative with respect to the spatial metric 7*^ . Note that Fourier transformation here is the 
generalization of Eq. Q for arbitrary spatial geometry of the universe. One can treat Q^j^x) as a plane wave in a 
fiat geometry case. 

Due to the existence of the second term on the left-hand side of Eq. (|D13p . we cannot solve this equation. Thus 
following |28| . we introduce the comoving momentum, q^^ = aP^. Regarding P as a function of comoving energy, 
q = g°, and conformal time, r, the third term in Eq. (|D5p may be replaced by ^§^ = ~hl^ij^^^''^ ^P ^^ ^^^ 
linear order. Then the linearized Vlasov equation, -^F{t, a;*, 7% g) = 0, becomes 

^/fc^i. f dFldhk ,„^„. 

_+zfc,.A, = ,_-_, (D16) 



where fk — fk (t, q, fJ.) ■ One finds the solution of Eq. IJD16|I as 

fk{T,q,ti) = e-^'''=(^-^'--Vfc(r.doc,<7,M) + §^ fdr' K{T')e-^^^^^^-^'\ (D17) 

where the prime on hki^r) denotes the derivative with respect to the conformal time. As there is no primordial tensor 
perturbations in the neutrino distribution function before neutrino decoupling, /fc(Tj,dcc, 9, m) = 0. 
The right-hand side of the linearized Einstein equation includes anisotropic stress as in Eq. (O ; 

<^-«' E /7|^n...g^,.(x), (DI8) 

where T^^ denotes the stress energy tensor of neutrinos. Since T,j^^ — —^=J—^q^qjF{q), its perturbation can be 
expressed as 

<57;M = «"' y ^ [*9J'5^ + i^l^^J + ^^5<ijW] ' (D19) 

/d\ 
— ^/,,fe(T,g,M)y7™QLW. 

The second and the third terms of l)D19p cancel out in linear perturbation theory. Thus 

n,,,Q^^.(x) = a-^ j'^qS'l'l'Tfx,kQL{^). (D20) 
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Inserting solution of the Vlasov equation (|D17|) into Eq. (|D20|l and using equality / dQqj'^j^j^j^e ''^'^'^Q 
1[8^SJ^ + (5*™(5J') / rf^ge-'^"(l - 2/i2 + /)g;^,„, one obtains 



= -4p-.(.) fdr' r#^^ ) /.L.(r'). 



r..AkHr-T'y 



(D21) 



Here, q^ = 097* and if — a qj^, and p^{t) 



'^ J d^qqF{q) is the unperturbed neutrino energy density, and a 



negative sign appears on the right-hand side of Eq. IJD21|) because integration by parts has been done. Also, we have 
used the identity 



^J dMl-V+/)e-^'"' = 



J2(w) 



(D22) 



Note that 



J2(-") _ J2 (") fOO h{u)j,, _ TT 
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^du 
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hju) 



_ , - , , _ J_ 3 

(-«)2 — ,j2 I J_oo «2 >*". — g , c^ii^i iiiil„^0 „2 — 15- 

Then the Einstein- Vlasov equation takes a form of an integro-differential equation; 



/i'fc(r) + fc'/ifc(r) = -24/,(r) 



a'ir) 
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Tv doc 
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(D23) 



and the fraction of the total energy density in neutrinos is 

Mr) 



Mr) 



-P{r) 



Vtyiao/a)^ 



MO) 



^Miao/a)^ + i^j + ri^)(ao/a)'^ 1 -t- aiT)/aEQ ' 



where 



MO) = 



n. 



i i-y I ^ ^ 1/ 



= 0.40523. 



(D24) 



(D25) 



The integro-differential equation (|D23|I was studied in j^^j 133. 13^4 l34l | in the cosmological context. Here we shall 
solve this equation numerically with all the Standard Model particles participating in the cosmic thermal plasma. 
Anisotropic stress, H^, vanishes during the matter era, as fi, —^ 0. Therefore, the damping effect is unimportant 
during the matter era. 
Following 1^, we write 



hxiu) = hxiO)xiu), 



(D26) 



which gives 



X"iu) 



2a'iu) 



x\u)+xiu) = -2AMu) 



a'iu) 



dU 

Uv doc 



J2iu-U) 
(u-[/)2 



x'(t^), 



(D27) 



where u = /ct, and derivatives are taken with respect to u. After the end of inflation, Tend, the amplitude of cosmological 
fluctuations is conserved until the mode re-enter the horizon, /ia(0) — /lA.kC'^'ond)- Note that the right hand side of 
Ea. (jD27|) disappears on the super horizon scales — neutrino free-streaming affects the tensor metric perturbation 
only inside the horizon. The initial conditions are taken to be 



X(0) = 1, x'(0) = 0. 



(D28) 



In the references |3 111331 . K[u\ 
i.e. K\u] = ^. 



sin u 3 cos 



zosu_ _|_ 3 sin It _ 1 (jq{u) + ^J2{u) + |j4{'u) j , which is the same function as our kernel, 
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FIG. 8: Comparison between numerical solutions and analytical solutions of tensor perturbations. The effect of neutrino 
free-streaming is included for numerical solutions, but not for analytical solutions. The dashed and short-dashed lines show 
numerical solutions of the high and low frequency modes, respectively. The higher fc-modes enter the horizon during the 
radiation era after neutrino decoupling, and thus the numerical solution is damped by neutrino free-streaming compared to the 
analytical solution, xi^''') = joikr) (solid line). On the other hand, the lower fc-modes enter the horizon much later, and thus 
the numerical solution is closer to the analytical solution during the matter era, xi^"^) ~ 3ji{kT)/kT (dotted line). 



We solve Eq. (|D27|I numerically by two steps; (i) we obtain a{T) and a'{T) from the Friedman equation (|26|l with (7*(t) 
in Sec. IIIII [Fig. 15]. and (ii) we solve Eq. (|D27|I with the scale factor that we obtained in the step (i) The numerical 
solutions as well as analytical solutions are presented and compared in Fig. |S1 The higher Fourier modes enter the 
horizon during the radiation era, but after neutrino decoupling. Thus they are damped due to the presence of the 
right-hand side of Eq. HU27|I . 

In order to estimate the damping effect, let us consider the radiation era after neutrino decoupling. During the 
radiation era, a'(u)/a = l/u, the analytical solution is given by x(m) = jo{u) in the absence of neutrino free-streaming 
in Eq. (|D27|I . In the presence of neutrino free-streaming, the solution becomes asymptotically (u 3> 1) 



X(u) 



A 



sin (u + 6) 



(D29) 



where A = 0.80313 and 5 = are obtained from our numerical calculations. This asymptotic solution is valid only 
for rather long wavelength modes which entered the horizon well after the neutrino decoupling. The suppression 
factor A^ = 0.64502 applies to the gravitational wave spectrum of the modes that entered the horizon after neutrino 
decoupling but before matter domination. 

In order to understand the shape of the spectrum. Figs. 01 and |S1 we need to consider shorter wavelength modes 
as well. This may be understood as follows. As we saw in Eq. (|D11|) , if the time derivative of the mode is negative 
(positive), the mode is damped (amplified). Integrating the amplitude of gravitational waves over time, the net effect 
of neutrino free-streaming almost always damps gravitational waves. This is because the contribution is mainly from 
the first period of x'(")j where the first trough is larger than the first peak. In previous paragraph we have considered 
the modes with fcr^dcc < 1- Now let us consider the higher /c-modes with fcr^doc '^ 1, or fc ~ 10~^° — 10~^ Hz. Note 
that fcr^dcc — 1 represents the mode which entered the horizon at the neutrino decoupling time, r^dcc- The mode with 
larger wavenumbers would enter the horizon earlier. Fig. [^ shows numerical solutions of x'('^) for which neutrinos 



decoupled at r^dcc given by kr^^e 



1.25, 2.5, or 3.75. For kr^^ 



1.25 and 2.5, neutrinos decoupled at the 



first trough of x'l"*^); where x'l"*^) is negative. Thus their amplitudes are damped by giving energy to free-streaming 
neutrinos (see Eq. (|D11|I and discussion below it). For fcr,ydcc = 3.75, neutrinos decoupled right after the first trough 
of x'(u), where x'('") is closer to zero. Thus its amplitude is unchanged, but its phase is delayed. Fig. 1101 shows 
numerical solutions of x'(") with kr^dec — 5.0. For kTi^^cc = 5.0, neutrinos decoupled at the first peak of x'(")j 
where x'i''^) is positive. Thus the amplitude of gravitational waves is actually amplified by gaining energy from free- 
streaming neutrinos, and we can see this feature on the spectrum. Fig. |31 at ~ 5 x 10^^*^ Hz. Neutrino free-streaming 
makes gravitational waves either damp or amplify depending on their frequencies. Note that this feature is generic to 
instantaneous decoupling of any kinds of particles, but not realistic for neutrinos as we mentioned in Sec. IIVI 
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FIG. 9: Derivatives of modes which entered the horizon before neutrino decouphng. The soUd Une shows an analytical solution, 
x' ~ —ji{u), during the radiation era without neutrino decoupling. The dotted, short-dashed, and dashed lines show numerical 
solutions of x'ikr) for which neutrinos decoupled at Ti/dec given by kr^dec = 1.25, 2.5, 3.75, respectively. They are damped by 
giving energy to free-streaming neutrinos. Vertical lines indicate the neutrino decoupling time for each mode. 



5 
^ 




FIG. 10: Derivative of a mode which entered the horizon before neutrino decoupling. The solid line shows an analytic 
solution, x' = ~ji{u), during the radiation era without neutrino decoupling. The dashed line shows numerical solutions of 
x'(fcr) for which neutrinos decoupled at Ti^doc given by fcr^dec = 5.0. The wave is amplified by gaining energy from free- 
streaming neutrinos. The vertical line indicates the neutrino decoupling time. 



For extremely short wavelength modes which have already been inside the horizon before neutrino decouphng, 
kr^dcc ^ 1 or A: > 10~^ Hz, the suppression becomes negligibly small; A ^ 1, but the phase delay, 5, is non-zero. 
These modes are undamped as positive and negative contributions of x' to the gravitational wave energy cancel out 
each other after several periods of x' ■ No net energy conversion from gravitational waves to free-streaming neutrinos 
would occur. 
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FIG. 11: The oscillatory factor [C + D] with respect to t = kn. The vertical line indicates s = kT2 — 300 < t and 
Ar = Ti — r2 = 0. The solid, dashed, short-dashed, and dotted lines show n = 0, 1, 2, and 5 respectively. The factor, [C + D]'^, 
takes on unity at Ar -^ regardless of n. 



APPENDIX E: OSCILLATION DUE TO DRASTIC CHANGE OF g,{T) 

In this Appendix we explain the effect on the gravitational wave spectrum from a sudden change in the number 
of relativistic species, g*. To do this, we need to calculate ri°i^*(A:2)/^"'(A:i), where ^2 ^ fci. In Sec. IIVI we have 
already seen the numerical prediction of the gravitational wave spectrum. In subsection IIII Bi wc provided the way to 
understand the relative suppression of gravitational waves at a given fc (= fci = ^2) with and without changes in g.^,. 
We shall discuss in a similar way what would happen to different Fourier modes, in order to fully understand imprints 
of (7* on the spectrum, such as oscillations, which are from cosmological events that change g■^, instantaneously or 
drastically. 

In Fig|31 we find an oscillatory feature at around 10"'' Hz, which corresponds to the mode entering the horizon 
at the QGP phase transition. At this energy scale, ~ 180 MeV, the effective number of relativistic species changes 
drastically, giving a sharp feature and oscillation in 0,^. To understand this, let us consider the simple analytical 
model employed in subsection IIII Bl Eq. 140|) is the mode which experienced such a change of g* and its coefficients 
A, B, C, and D are 



A{s,n) 



(El) 



^-^ [ - 2sYi+yiqr4;^/2(s) sins + F/T+4;r/2(s) (-2scoss + (3 + Vl + 4n) sins) ] , 
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4s3/2 



2s-^i+VT+4;r/2('S) sins + J,/t+^,/2{s) (-2scoss + (3 + Vl + 4n) sins) 



(E2) 



C{s,t, n) 



A^/7t 



;n7r[- 2iJ_„_3/2(i) cosi (j„+i/2(s)(s coss - sins) + sJ„+3/2(s) sins) 



-2J_„_i/2(i) (J„+i/2(s)(scoss- sins) + sJ„+3/2(s) sins) (cost + isini) 
+2 (j_„_i/2(s)(scoss — sins) — sJ_„_3/2(s) sins) 
(-t J„+3/2 (t) cos i + J„+i/2 (i) (cos i + i sin t)) ] , 



(E3) 



D{s,t,n) = 
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Wst 



secn7r[2tJ„_(.i/2(0 cosi (J_„_i/2(s)(— scoss + sins) + sJ-„-3/2(s) sins) 



-2J_„_i/2(i) (J„+i/2(s)(scoss- sins) + sJ„+3/2(s) sins) (tcost ~ sini) 

-2sini(iJ_„_3/2(t) (j„+i/2(s)(s cos s - sins) + sJ„+3/2(s) sins) 

+ {-tJn+3/2(t) cost + Jn+i/2{t)) (s J_„_3/2 (s) sin s + J_„_i/2 (s) (-S COS s + sin s)) )] , 



(E4) 



where s = fcr2, i = fcri and s < t. Here, Jn{x) and yn(a;) are the Bessel function and Neumann function, respectively. 
At this time, we are interested in different fc-modes, fci < /c2- (However we evaluate ^h{k) at the same time, r.) We 
obtain 






fc|/l2(fc2T) 

kjh^{kiT)'' 
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C(fc2) 



JoihT) 



D{k2) 



yoihr) 



jaiki-r) 



2n r 



C(fc2)|+i?(fc2)|^ 



2n 



= /£ [c(fc2)+i?(fe)]^ 



(E5) 
(E6) 



where ~ denotes the subhorizon limit, « denotes the asymptotic limit as kr —> large, and — > denotes the limit in 
Afc = fc2 — fci -^0. Eq. IIE5|l tells us the exact ratio between different /c-modes. While we obtained only the 
suppression factor, {t2/ti) ", in subsection IIH Bl we now also obtain the oscillatory factor, [C + D]"^. Fig. 1111 shows 
that the factor, [C + D]"^, oscillates and takes on unity at At -^ regardless of n. Here, n = 5 represents g* ex r^'^", 
which is an extremely drastic change. This gives us a complete analytical account of the shape of Fig. |5| 
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